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Abstract 

In their textbook, Suzuki and Varga [Y. Suzuki and K. Varga, Stochastic Variational Approach to Quantum- Mechanical Few- 
Body Problems (Springer, Berlin, 1998)] present the stochastic variational method in a very exhaustive way. In this framework, 
the so-called correlated gaussian bases are often employed. General formulae for the matrix elements of various operators can 
be found in the textbook. However the Fourier transform of correlated gaussians and their application to the management of a 
relativistic kinetic energy operator are missing and cannot be found in the literature. In this paper we present these interesting 
formulae. We give also a derivation for new formulations concerning central potentials; the corresponding formulae are more 
efficient numerically than those presented in the textbook. 



PACS numbers: 02.70.-c, 12.39.Pn 

I. INTRODUCTION 

There exist several different technical methods to solve 
the few-body problem with accuracy: Monte Carlo cal- 
culations Faddeev and Yakubovsky treatments [2(], 
hyperspherical formalism Q , expansion on various types 
of orthogonal Q , or non orthogonal bases [f| . Each tech- 
nics shows specific advantages and drawbacks. Among 
others, the stochastic variational method is especially at- 
tractive. It relies on expansion of the wave function in 
term of gaussian type functions. The stochastic algo- 
rithm allows to consider very large bases with a minimum 
of variational effort. The drawback of this method is the 
non orthogonality of basis wave functions with the pos- 
sibility of appearance of spurious states due to overcom- 
pleteness; if this last inconvenient is overcome, using non 
orthogonal bases is not really a problem. The generalized 
eigenvalue problem arising in this case is well under con- 
trol nowadays. The great advantage of using gaussian 
type functions is the rapid convergence and, above all, 
the possibility to compute the resulting matrix elements 
with analytical expressions most of time. 

The stochastic variational method is described in full 
details in the remarkable textbook by Y. Suzuki and K. 
Varga Q , where most important and fundamental formu- 
lae are derived. This very complete work will be refered 
as SV throughout this paper, and all subsequent refer- 
ences can be found in it. 

However, working on other projects, we were faced to 
the necessity to use some particularly important matrix 
elements that arc not found in the exhaustive SV text- 
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book. In particular, the Fourier transform of the corre- 
lated gaussians and the matrix elements of a relativistic 
kinetic energy operator are dramatically missing. In this 
paper, we want to complete SV with important formu- 
lae that are of crucial importance for some applications. 
Besides the relativistic kinetic energy, we derive also new 
expressions for the matrix elements of central potentials. 
These expressions are simpler and more efficient numer- 
ically than those given in SV. We also report them here. 

To achieve some unity, to precise our notations and 
to have a self-contained paper, we will also present be- 
low some formulae already present in SV. In this case we 
give the references where they can be found in this work. 
We will derive expressions for the most general correlated 
gaussians (arbitrary number of particles AT + 1, arbitrary 
angular momentum L) but we restrict ourselves here to 
operators that do not mix spatial and spin degrees of 
freedom, so that we are concerned only with the ma- 
trix elements for the spatial part of the wave functions. 
Moreover, the non-natural parity states are very difficult 
to handle in correlated bases (this is possible but the ex- 
pressions are much more involved) and, in the following, 
we just study natural parity (i.e spatial parity equal to 
(-1) L ) states. 

The paper is organized as follows. We first describe 
the systems under consideration (intrinsic coordinates, 
definition of correlated gaussians and their generating 
functions). Then, we give the matrix elements for the 
overlap, the non-relativistic kinetic energy and a first ex- 
pression for central potentials. Except for special points, 
most of these formulae can be found in SV. Lastly, the 
novelties concern new expressions for central potentials, 
the Fourier transform of correlated gaussians and its ap- 
plication to relativistic kinetic energy. 
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II. THE SYSTEM UNDER CONSIDERATION 

The stochastic variational method can be applied to 
systems composed of more than one particle up to around 
ten particles. Let us denote by N + 1 the number of 
particles (N > 1). The particle i is located at Pi, so that 
the position vector relative to some origin O is OPi = Ti, 
and has a corresponding conjugate momentum pi. Since 
we are concerned only with spatial degrees of freedom, 
we ignore the color, isospin and spin variables. 

It is interesting to introduce the center of mass coordi- 
nate Rem = xn+i and the corresponding total momen- 
tum P = tvn+i- The intrinsic description is expressed in 
term of internal coordinates. There exist different pos- 
sible choices depending upon the nature of the studied 
systems. Here, we are concerned only with the Jacobi 
coordinates. 



A. Jacobi coordinates 

In order to simplify the notations, let us note mx2,„% = 
nil + m2 + . . . + wij and Gi the center of mass of the first 
i particles: OGi = (m^i + m 2 r 2 + . . . + m i r l )/mi2... l . 
The Jacobi coordinate x. L is defined as the position of 
the particle i + 1 relative to the center of mass Gi of the 
previous particles. Explicitly 



It is a matter of simple calculation to check that the 
matrix elements of the inverse matrix are given by: 



T _ u __ _ mj±L _ k<l< v 



mi2...i+i 

TOl2...j 



I < N, 



(u- 1 )^ = k>l + l, 

(U-^kN + x = 1 Vfc. 



(•5) 



The conjugate momentum for variable Xi is denoted 
TTi = —id/dxi (h = 1) and must fulfill the conditions 
[(xi)k, (iTj)i] = iSijSki- It is easy to check that the rela- 
tions corresponding to Eqs. write 



JV+l 

3=1 
N+l 

i> X! '■•' 7r .'- 

3=1 



(6) 
(7) 



With this definition, the total momentum P = pi + p2 ■ 

. . . +PN+1 IS just 7TJV+1. 



B. Correlated gaussians 



Xi = OGi-r i+1 , (1) 

There are obviously N Jacobi coordinates Xi (i < N) cor- 
responding to intrinsic coordinates and the special vector 
xn + i corresponds, as already mentioned, to the center of 
mass coordinate for the system (this vector has a physical 
meaning for particles with non vanishing masses; for the 
case of null mass particles, a way to get rid of this prob- 
lem is to employ a relativistic kinetic energy, as explained 
in section IVIIII) . 

It is important to express the Jacobi coordinates in 
terms of the original position vectors 

N+l 

Xi=J2 U n r r ( 2 ) 

3=1 

It is easy to calculate, from definition (TTJ), the value of 
the elements Ua : 



Wi2...i 

Uu+i = -1, (3) 
Uij =0 if j > i + 

The inversion of relation ([2]) is often useful 

N+l 
3=1 



In the stochastic variational method, the wave function 
for the system is expanded on gaussian type functions. 
There exist several different types of such basis states. 
Here we are concerned only with the so-called "correlated 
gaussians" . The space function is expressed in terms of 
the Jacobi coordinates. 

The simplest version is a function that is a product of 
separate gaussian function for each Jacobi coordinate. In 
this case the calculations are very easy, but the conver- 
gence is slow, and, moreover, the available Hilbert space 
is limited, because each pair of particles is in a S-state. A 
more elaborated version considers instead the argument 
of the exponential as a bilinear combination of the Jacobi 
coordinates: Ylij=i ^ij x i • Xj. The matrix A must be 
symmetric (A — A) and positive definite. Obviously, the 
Hilbert space is enlarged and the convergence is acceler- 
ated. 

To deal with a non vanishing total angular momen- 
tum, one must introduce spherical harmonics somehow 
or other. The most elegant manner is to use a single 
solid harmonic D^iiw( u ) = V l Ylm(v)- To achieve some 
symmetry, and also to have more variational parameters 
at our disposal, the argument of the solid harmonic is 
the most general linear combination of the Jacobi coor- 
dinates V = J^iLl u i x i- 

In order to simplify the notations, let us introduce a 
"super- vector" x = (xx, X2, ■ ■ ■ , a; at) and write the co- 
efficients of linear combinations as a "line (or column) 
matrix", i.e u = (ui, U2, • • ■ , u n ). This allows to shorten 
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the expressions using the usual matrix operations. For 
example (the presence of the symbol • deals with a spa- 
tial scalar product, while the absence deals only with a 
linear combination) 



N 



N 

E 



(8) 



With those definitions, the most general correlated 
gaussian is given by (note a slight difference with SV 
notations; their matrix A is twice ours and moreover TV 
is the number of Jacobi coordinates while SV consider it 
as the number of particules) 

fKLM{u,A;x) = exp(-x ■ Ax) \ux\ 2K y LM (ux). (9) 

The correlated gaussian represents the basis state in coor- 
dinate representation fx lm{u, A; x) = {x\ipKLM (% A)). 
Thus each basis state is described by N(N + 3)/2 free 
parameters (N(N + l)/2 for the matrix A and N for the 
vector u). This prescription (J9j) is only able to deal with 
natural parity states. The term |ua;| 2if is introduced for 
generality and to treat with more accuracy potentials 
with specific singular features. However, it complicates 
a lot the resulting expressions. It is often more conve- 
nient (except when the potential is so singular that the 
resulting integrals diverge) to keep in the calculation the 
correlated gaussians restricted to K — 0, including more 
basis states to compensate a slower convergence. The 
even simplest version corresponds to a diagonal A ma- 
trix. 



the correlated gaussians since one has 
fKLM(u,A;x) 



B 



KL 



( d 2K+ 

x / de Y LM (e) I gx2K+L g(^eu, A: x) 



l , (11) 

A=0,|e|=l 

where the geometrical coefficient Bkl is defined as 

4n(2K + Ly. 



B 



KL 



2 K K\ (2K + 2L + 1) 



(12) 



In Eq. (jlip . the super- vector s = Xeu must be under- 
stood with all its components proportional to the same 
three vector e, namely = Ai^e. 

Using Eq. (JTTJ) in the expression of the searched matrix 
element leads to 



(ipK'L'M'iu ', A')\0\lpKLM(u, A)) 

' r dede'Y LM {e)Yl, M ,{e') 



Br'L'Bkl 

q2K'+L'+2K+L 
dx ,2K>+U dx2K+L 



(13) 



(O) 



/ A=A'=0,|e| = |e'| = l 

with the matrix element between the generating function 

(O) = <s(AW, A': x)\6\g(Xeu, A; x)) (14) 

The matrix element that is left for computation is now 
between the generating functions, the form of which is 
much simpler. 

In the rest of the paper, we focus our study on oper- 
ators that are scalar for spatial coordinates, so that the 
only non vanishing elements are those with L = L', and 
M = M'. Moreover, the matrix elements do not depend 
on the magnetic quantum number M. 



C. Matrix elements and generating functions 

We are interested in the calculation of the matrix ele- 
ments for some operator O on the correlated gaussians, 
namely {ipx'L'M 1 (u' , A')\0\iPklm(u, A)). The critical 
point in the computation of such element is to manage 
the difficulties due to the presence of the solid harmonics. 
An artful way to deal with them is to use the generating 
functions for the correlated gaussians. Since the techni- 
cal details can be found in SV, we just recall below the 
most important results and let the reader have a look on 
SV (in section 6.3) for rigourous proofs. 

Let us define the functions 

g(s, A; x) = exp(— x ■ Ax + s ■ x) (10) 

where s is an arbitrary super- vector, s = («i , S2, • ■ ■ , sjsr), 
and, as usual s ■ x = J2i= 

The g functions are called the generating functions for 



III. OVERLAP 

The overlap between basis states Mk'kl — 
(ipK'LM(u',A')\ip KLM (u,A)) is a crucial ingredient 
in the equation of motion. Since the basis states are non 
orthogonal, there is no reason that such an element is 
diagonal. 

As already mentioned, it is easier to calculate first the 
elements for the generating functions. The technics is 
always the same; we first diagonalize the matrix A and 
change the integration variables to the eigenvectors of A. 
More details are presented in SV, and we will develop 
more deeply this technics later on. Explicitly, one finds 
(see SV, Table 7.1 page 124): 

/ N \ 3 / 2 

(g(s',A';x)\g(s,A;x)) = f^g) ex P(^ ' B~ 1 v) 
= Mo (15) 
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where 



A simple calculation leads to the explicit expression 



B = A + A'; v = s 



(16) 



The overlap is obtained introducing result (|15l) in the 
general formula (fHij). Since the variables A,A',e, e' ap- 
pear only in v, the term depending on det B factorizes 
out the integral. The exponential is then expanded as 
a series in A, A',e • e'. The derivation and integration 
are performed without difficulty. The final result is (SV, 
(A.6) page 248) 



TnR = T — T cm 
N 



where the super-vector -k = (tti, 7t 2 , ■ ■ . , ttn) is intro- 
duced. The matrix A is symmetric and its general ele- 
ment writes 



M K , KL ^ K 'l L ^ K + L ^ (^T (17) 



Br'lBkl 



detS 



min(K,K') 



k=0 



n K—k iK'-k 2fc+L 

(K — k)l (K'-k)\ {2k + L)\ 



with the following complementary definitions 

q = ^uB- l u] q' = -u'B-V; 

p = -u'B^u = -uB^u'. 
H 2 2 



(18) 



For the important peculiar case K' = K = 0, this formula 
simplifies a lot and we are left with (SV (A. 7) page 249): 



4?r \detB 



N+l 



Ay — ^ UikUjk ■ 

7 -, ^k 



(21) 



fe=i 



The matrix elements of Tmr on correlated gaussians 
rely again on the expression of the elements on the gen- 
erating functions. This calculation is a bit more com- 
plicated, but one can show that (see SV, Table 7.1 page 
124): 



{g(s',A';x)\7t ■ A-K\g(s,A;x)) 



Mo [6Tr(AB _1 A'A) - y ■ Ay] 



(22) 



where Mo is precisely the overlap expression (|15p . Tr(A) 
means the trace of the X matrix, and the super- vector y 
is given by 



y = A'B^s - AB^s'. 



(23) 



IV. NON RELATIVISTIC KINETIC ENERGY 

It is well known that the use of Jacobi coordinates 
allows to share the total kinetic energy operator T — 
J^iL^ 1 Pi/(^ m i) m t° the kinetic energy of the system as 
a bulk T cm — P 2 / (2M) (M = mi2...N+i) and an intrinsic 
kinetic energy operator T^r depending only on intrinsic 
momenta. 



It remains to insert result (|2"2"|) into Eq. (fT3"|) . The 

variables A, A', e • e' appear now still in the exponential 
part of Mo but also in y ■ Ay while the trace is indepen- 
dent of these. Fortunately, the new contribution is just a 
polynomial in term of these variables. Thus, developping 
the exponential as a series leads to a total contribution 
which is polynomial and which can be treated exactly in 
the same way as the overlap. 

The final result is (see SV, (A. 10) page 250) 



(4>k'lm(u , A )\tt ■ Att\iPklm(u, A)) 



(2K' + L)\ {2K + L)\ ( ir N x 3/2 



Bk>lBkl \AetB t 

™in(K,K') K-k-1 iK'-k-l 2fc+L-l 

x J2 BkL [Rqqp + P(K - k)q' p + P'(K' - k)qp + Q{2k + L)qq] ^ _ fc) , H — - — ^ - — (24) 



with the numbers P, P' ,Q, R defined by Again, the formula for the special case K = K' = is 

P = -uB- 1 A' AA 1 B'^u ; P' = -u> B^ 1 AAAB^ 1 u '; 
Q = 2u'B~ 1 AAA'B~ 1 u ; R = 6Tr(AS _1 A'A). (25) 
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much simpler 

V p 

(26) 

The pleasant feature is the factorization of the overlap 
matrix. 

Clearly the matrix U, transpose of U defined in Eq. 
©, is different from U^ 1 defined in Eq. ©. Thus U 
is not an orthogonal matrix. Nevertheless, one has the 
property detU = 1. This matrix has also another very 
interesting property. 

Let us define the reduced mass fii for particle i by 



f 

(H 



1 



1 



77112. 



rrii+i 



(27) 



It can be shown that the definition ([2|) allows to write the 
non-relativistic kinetic energy operator in a form without 
mixed terms, namely 



1 



N 9 
x - 7T,- 



T NR = -fr-A7r = ^^. 



(28) 



This means that, in the matrix elements (j2"4")) - (f2"6")) . one 
can restrict the calculations to a diagonal matrix A since, 
in this case, Ay- = (l/fii)8ij. 



V. CENTRAL POTENTIALS: A FIRST EXPRES- 
SION 

In the most interesting cases, the potentials appearing 
in the few-body problem are either one-body potentials 
or two-body potentials. In the first situation, the poten- 
tial has a form V = ^ Vi (r< — R cm ) while the second 
situation the potential has a form V = J2i<j ( r i ~ r i)- 
This is the consequence of translational invariance. In 
both cases, the argument of any potential is always a 
linear combination of Jacobi coordinates, so that the 
most general form of the potential is a sum of terms like 
V(wx). 

In this case, one can write generally 

(i/>\V(wx)\tl>) = J V(r)(ij\5 3 (wx-r)\^) dr. (29) 

It is necessary to calculate the matrix elements of 5(wx — 
r) on the generating functions; these are given in SV 
Table 7.1 page 124. 

However, we are interested here in central potentials 
that are also invariant under rotations; the argument of 
the potential must be now \wx\. In this case, the equiv- 
alent expression of (|29|) needs a single integral 



(i>\V(\wx\)\ip) = J V(r) (ip\ S(\wx\ - r) dr. (30) 

It is necessary to calculate the matrix elements of 
<5(|u)a:| — r) on the generating functions. They are not 



given in SV. It is possible to get them either directly by 
elementary integrations, or by the developing 6 3 (wx — r) 
in spherical coordinates and integrating on angular vari- 
ables. We find 



{g(s',A';x)\5(\wx\-r)\g{s,A;x)) = 



(wB^wf/ 2 



'o 



r\wB v\ 
wB^w 



exp 



AM 



r 2 + {wB^vy/A 



wB~ 1 w 



(31) 



where io(z) — sinh(z)/z is a modified spherical Bessel 
function. 

One must implement the element (|31[) into Eq. JT3]). 
For further convenience, we give here more informa- 
tion concerning this computation. Let us note c = 
2/(wB~ lr w). First, it is possible to remove the term 
exp(— cr 2 /2) and all terms depending on c only out of 
the derivation and, even, out of the integral. The terms 
containing the variables A, A', e • e' are present in the ex- 
ponential part of M.q which writes explicitly 



M oc exp(<jA 2 + q'\' 2 + pXX'e ■ e'). 



(32) 



These terms also appears in the variable z = (wB~ 1 v)/2. 
Let us introduce the numbers 7 = c(wB~ 1 u)/2, 7' = 
c(wB~ 1 u')/2. The variable z now reduces to cz = 7Ae + 
7'A'e'. Thus, the term zo(cr|z|) exp(— cz 2 /2) must also 
be included in the derivation term. 

Using the fact that exp(— s 2 + 2sx) is the generating 
function for Hermite polynomials H n (x) , it is possible to 
prove the interesting relation 

{Az 2 ) n 



io(Bz)eM-Az') = 



B ^ (2n + l) 

n=0 V ; 



(J ff 2 „ +1 (73/2VA). 

The technics presented in SV consists in employing 
this formula to evaluate io(cr|z|) exp(— cz 2 /2). Forget- 
ting terms depending on c only, there appears a term 
depending on r: H 2n +i(r y/c~/2) / r and a term 



(c 2 z 2 Y = ( 7 2 A 2 + 7 ' 2 A' 2 + 277'AA'e • e') r 



(34) 



which must be maintained in the derivation and which is 
expanded as a series. It is gathered with the series coming 
from exponential |32|) . After this grouping, the general 

monome of this series has the form A S A' S (e • e')*. The 
derivation Q 2K + L /d\ 2K+L being taken for A = leads 
to a term (2K + L)\ S2K+l, s ', one has also a term {2K' + 
L)\ S2k>+l,s> coming from the derivation with respect to 
A'. Lastly, the term (e-e') 4 is expanded as (do not forget 
that |e| = 1 = |e'|) 



(e-ey^TT^T — 



t\ 



{t-l)\\{t + l + 



Y; m (e)Y lm (e') (35) 



and the integration over e and e' leads to a term 

BkL $2k+L,t- 



5 



After this tedious but straightforward cooking, and us- 
ing expression ([30|) . we get the final result (SV, (A. 128) 
page 282): 

(lljK'LM{u',A')\V(\u!x\)\tp K LM(u,A)) = 

ir N \ 3/2 {2K 1 + L)\ (2K + L)_! A '1g +L c ) 



detB/ B K , L B KL ^ 

min(K,i<r') 

x ^ B kL F^_ ktK/ _ k!2k+L (q,q' >P> -/,Y) 

k=0 



(36) 



In this formula, the dynamical quantities q,q' , p have be 
defined in Eq. (|18p for the overlap, while the new ones 
c, 7, 7' are specific to the potential 



wB 1 w 1 



7 



wB l u , w_B 1 m 



wB 1 w 



; 7 



— (37) 



The form of the potential appears explicitly through the 
following integral 



J(n, c) 



Lastly, the geometrical function F is defined as 



FZ p , il (q,q',p, % i) = n\ £ 



(38) 



(39) 



qP -m q iP ~m 
(p — to)! (p' — to')! 

pl — n-\- r m-\-rn' ryn+m — rn'^l n ~ ra + rm 

(I — n + to + to')! 2 m+m ' ml m'!(n — to — to')! 



We tried to express this function in term of hypergeo- 
metric functions, but without success (when no explicit 
bounds arc indicated, the indices in the summations run 



on values that do not give integer negative values for the 
factorials). 

The expression (|36]l has the nice property to allow 
easy checks. The bra-ket symmetry is trivial. For such 
a symmetry, K <-> K' , A <-> A', u <-> u', q <-> q', 
7 <-> 7', while -B, p and c remain unchanged. The sym- 
metry follows from the property F™ , ; (g, g', p, 7, 7') = 
Fp, pl (q',q,p,y,-f). It is also easy to recover the over- 
lap expression if we set V(r) = 1. In this case, due to 
the orthogonality of Hermite polynomials, one has n = 0, 
J(n, c)/c" = 1, F p ° )pV = qPq' p ' p l /[p\p'l 11} and the overlap 
expression follows immediately. 

As always, there is a great simplification in the special 
case K = K' = 0. Indeed, one has the simple expression 
(SV, (A. 130) page 282) 

(^lm\V(\wx\)\^ lm) =M L L\J2 7T^4t (—) ■ 

^ {L- n)\ \pc J 

Here again, the overlap factorizes. 

A small drawback in the formulae (|3"6")) and is the 
presence of the Hcrmitc polynomial H^n+iix). One way 
to deal with it is the use of recursion formulae. How- 
ever, the accuracy decreases with increasing n; more- 
over, closed expressions for J(n, c) do not exist for various 
forms of potential V(r). There is a way to get rid of these 
two difficulties. 

We define new integrals: 



V{u)u k e 



-Au* 



du. 



(41) 



Indeed, a lot of closed expressions exist for various forms 
of potentials V(u). Expanding the Hermite polynomials 
as a series, and rearranging the summations leads to a 
new expression for the matrix elements in terms of Ty 
functions. This form is absent in SV. 



&k'lm(u ,A r )\V{\wx\)\ip KLM {u,A)) = l \ 



2ir2 K + K '- 1 B K , L B KL \detBj (- C ) K + K '+ L 

H--H--L . „ ..„ ,-- ..„ , n _ /n \ m HK,K') . ( 42 ) 

X 

n=0 ' fc=0 



+f +L (-2c)"^ y (2n + 2,c/2) ) KjK ,, L (2rf V7 pc\ 

h ( 2 « + !) ! h V PI ' rf ' 77' ) 



with function H defined by and 



s\(K - k - sY.s'HK' - k - s')\ 

V ( if ^ + X +L-r)ly K,K',L ( ,s 1 

h ) (K + K' + L-n-r)^ [X ' X) * (r - s - s')\{2k + L + s + s' - r)\ ' 

(43) (44) 
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Due to the property G%'£ r ' L (x,x') = G^^'' L (x',x) the 
hermiticity of the element is transparent. 

For the special case K = K' = 0, this formula reduces 

to 



yK',K,L, i 



{^ 0LM (u',A')\V(\wx\)\^ QLM (u,A)) =Af L 2c x — LI 

Ztt 



xV^t,^ r^Vfl-^V'". (45) 



H (2n + !)!(£ ■ 



This formula needs essentially the same numerical effort 
than expression (|40p but it is written in terms of a much 
more convenient integral. 



VI. CENTRAL POTENTIALS: A NEW EXPRES- 
SION 

In this section we present new formulae for the ma- 
trix elements of a central potential; the corresponding 
expressions are more efficient for a numerical treatment. 

Let us remark that the relationship (f3~4")l gives a link 
between e-e' and the variable z. This opportunity allows 
to put the exponential (|32[) under the form 



exp(qA 2 



q X 



pc 2 z 2 /{2 11 1 )) 



(46) 



with the new variables q = q — pj/(2j') and 
q' = q' - pi/{2 1 ). 

One gathers the part of exponential depending upon z 
with the other exponential depending upon z into a single 
exponential which takes the form exp(— cz 2 /2) with the 
definition c = c(l — pc/^j')). We then treat the term 
io(cr|z|) exp(— cz 2 /2) with the already mentioned rela- 
tion ([33]) . The rest of the derivation is essentially similar 
to that of section [V] The trick of gathering two expo- 
nentials into a single one allows to gain one expansion 
into a series; the price to pay is the use of renormalized 
dynamical quantities. 

The final formula, which is absent in SV, looks like 

(iPk>lm(u',A')\V(\wx\)\^ klm {u,A)) = 



(2K' + L)\ {2K + L)\ ( 

Bk'lBkl \detB 

K+K'+L 

x E (yc) j K a > c ) 

mm{K,K') 2k+L 



(47) 



n t?K,K'.L, / / \ 



k=0 



In this formula, the dynamical quantities c, 7, 7' are iden- 
tical to those defined in Eq. (|37f while new quantities are 
necessary 



q = q 



PI 
27" 



PI 
2 7 



1 



PC 

77' 



(48) 



The form of the potential appears explicitly through an 
integral of new type 



J(n,a, c) 



1 



V^(2n + 1) 



(49) 



x / V(xy/2^)e~ ax H 1 (x)H 2n +i(x) dx. 



Lastly, the geometrical function F is defined by 



F n k {x,x ,y,y) = n\ 



min(«-fc.A'+L) ^K+L — m 



E 



m — max ( k + L , n — K ' ) 



(K — n+rn 



y2m-L^/2fn-m)+L 



(if — n + m)! (m — k — L)\(n — k — m)! 



(50) 



Comparison between equivalent forms (|36[) . given in SV, 
and (|47p . novelty of this paper, warrants some comments. 



• They both have very similar aspect. The new form 
requires the calculations of 3 new quantities q,q' ,a 
and an integral J(n, a, c) instead of J(n, c). But in 
doing this, the numerical effort is essentially the 
same. In contrast, the rest of the calculation is 
much more efficient under the new version; indeed 
the F function depends on 4 variables instead of 5, 
it requires a single summation instead of a double 
and the general monome needs one factorial less. 

• The equivalence between both can be proved using 
the following very interesting relationship 



, 2 1 Tvl / f(y)e-^H 1 (y)H 2L+1 (ay)dy 



L , 2 



E 



(2L+1)! 



a 2 -I) (£-n)!(2n + l)! 



(51) 



x / f(y)e-y-H 1 (y)H 2n+1 (y)dy 



valid for any function f(y). 

• The symmetry properties are more transparent in 
the SV version, but they can be proved as well 
in the new version. Once the equivalence between 
both is shown this point is of minor importance. 

The new version is especially interesting in the peculiar 
case K = K ' = since the result writes 



(^olm\V(\wx\)\^ olm ) =N l L\ 



l-a 



In this case, there is no summation at all! 



a 3/2 J(L,a,c). 

(52) 
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To finish this part, let us derive the general formula 
expressed in term of the T integral instead of the J in- 
tegral. The calculation relies essentially on the technics 
presented in section [V] The result is 



(■>Pk'lm(u', A')V(\wx\)\ip KLM (u, A)) 



4{2K' + L)\ (2K + L)\ ( en 



N x 3/2 



V^Bk'lBkl \2detB 

K+K'+L 



E 

n=0 



F v (2n + 2,c/2) 



(53) 



(2n + l)\ 



min{K,K') 2k+L 



V — a..nK,K ',L,- -, „, 



:B kL G k ' (g,g',7,7',c) 



with function G defined by 

GK.K' .L / - -i / \ 



(54) 

^ gJC+L-r^/AT'+r-n 7 2r-Ly2(n-r)+L q^ 2 

(tf + Z,-r)!(r- ifc-L)! n ' fe ' r( 2cJ ' ^ 



(s + rc)! 



and with 

A w +r — n 

M «,fc,r(^) =E (K' + r- s - „ i!(.s f „ /, /■)! .,•! ' 

s=max(0,fc+r-rf) ' v ' 

This formulation is also more efficient numerically than 
(142p for the same reasons as those explained above. 

Application of this formula to the special case K — 
K' = leads exactly to the formula lf45|) got before and 
thus does not bring any novelty. 



VII. FOURIER TRANSFORM OF CORRE- 
LATED GAUSSIANS 

For a number of applications, it is useful to have the 
Fourier transform (FT) of correlated gaussians. This 
point is missing in SV and, to our knowledge, such an 
expression seems absent in the literature. In this section, 
we want to derive the corresponding relation. 

The basic ingredient for this calculation is the FT of 
a correlated gaussian limited to one variable, namely the 
value of the integral 



(27T) 



-3/2 



-iqr 



e ar r 2K 



yLM{r) dr. 



(56) 



The first thing to do is to use the traditional devel- 
opment of a plane wave in terms of spherical harmonics 



e - iq .r = 4 tt S£o (-iYji (qr) S+L-, Y lm (g)F^ (f) where 
ji(z) = y / Tr/(2z)J l+1 / 2 (z) is the spherical Bessel func- 
tion. Integration over angular r variables leads to a 
&im,LM factor which restricts the infinite summation to 
the single term I = L,m = M. Apart from some 



constant factors that can be gathered outside the inte- 
gral, we are left with a radial integral, which, after a 
trivial change of variable, can be put under the form 
/ °° e -x 2 x 2K + L +z/ 2 J L+1 / 2 {qx/^/a). Fortunately, this in- 
tegral is analytical and can be expressed in terms of gen- 
eralized Laguerre polynomial L^. This miraculous for- 
mula can be found in Ref. @ (formula (6.631) n 10, page 
738); it writes (it can be obtained using the generating 
function of Laguerre polynomials expressed in term of 
Bessel functions 0) 



e- x2 x 2n+t2+1 J^(2x^)dx = ^ e - z z^ 2 L»(z). (57) 



In this formula fi is any real number > — 1 and n is a 
positive or null integer. The case n = is simpler due 
to the property Lq(z) = l,Wz. We thus get the searched 
integral 



1 



(2tt) 3 / 2 



-iqr -ar 2K 



r 2K y LM (r)dr 



«_I_ l l+i/2 / r_ 



{-if 



4a J (2a) L + 3 / 2 



V/(4Q) 3W<?). 



(58) 



To get the formula corresponding to K = 0, it is sufficient 
to replace in Eq. (|5"8|) the terms depending on if by 1, 
at the head of the right hand side expression. 

To derive the general FT, we begin, for pedagogical 
reasons, with the special case K = 0. Let us note 



h 0LM (u,A;ir) = (2tt) 



-3N/2 i e ~ l7r ' x f LAi{u, A; x) dx 



= (2^r 3jy / 2 



X-Ax — ITT-X 



y L M(ux) dx = I. (59) 



The matrix A is symmetric and thus can be diagonal- 
ized with an orthogonal matrix T : A — TDT. A being 
positive definite, all the eigenvalues of the diagonal ma- 
trix D are positive, and there is no problem to build the 
square root D 1 / 2 of this matrix. Instead of the original 
x variables, it is wily to work with the new variables 
z = D 1 / 2 Tx. At this stage, one has 



I = (27T) 



-3AT/2 



(det A) 



-3/2 



dz exp(-i • z — zir • TD- 1/2 z)y LM (uTD- 1/2 z). 

(60) 

In a second step, we change again from z to new Z vari- 
ables with the help of an orthogonal matrix U : Z = Uz. 
The Jacobian of the transformation is unity, and, more- 
over, z ■ z = Z ■ Z. There is a large freedom in the choice 
of the U matrix. One can take this opportunity making 
a choice such that the argument of the solid harmonic is 
proportional to one of the new variables, let us say Z\. 
Now, the integral takes the form: 



7=(2^)- 3Ar / 2 a- i (detA)- 3 / 2 



x / dZ exp(-Z ■ Z - iV ■ z)y LM (z 1 ). 



(61) 
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with a = {uA^u)- 1 ' 2 and V = UD- l l 2 fn. 

Since the integral is separable for each variable, the rest 
of the proof is quite easy. N — 1 integrals, corresponding 
to Zk, k > 1 are equal to 

(2ir)-V 2 J dZ k cxp(-Z 2 -iV k -Z k ) = 2- 3 / 2 exp(-V fe 2 /4) 

(62) 

while the integral relative to Z\ can be calculated using 
expression ([58)) to get 



(2tt)- 3 / 2 J dZ l exp(-Z 2 - iVt ■ Z x )y LM (Z x ) = 

{ _ i} l 2 -l-3/2 exp (-y 1 2 /4)3; iM (y 1 ). (63) 

Grouping all the exponentials into a single one gives an 
exponential with the argument — (V 2 + V 2 + . . . V^)/4 
= -(V ■ V)/A = -(tt • A- 1 tt)/4. Thus, the final result 
for the FT looks like 



h aLM {u,A;ir) = 



(-i) L e~i*- A 
2 L + 3N / 2 (det Af/ 2 



y L Ai{uA 1 tt). 



(64) 

It is remarkable that, in this special case K = 0, 
the FT of a correlated gaussian is proportional to 
a correlated gaussian since we have the property 

{-i) L , A^ 1 

h 0L M(u,A;ir) = 2L+3N/2 ^ detA ^ /2 foLM(A^ 1 u,—;7r). 

(65) 

With this relationship, it is instructive to check the con- 
servation of the norm, namely Ml = (JolmI/olai) = 

(hoLM\hoLAl) ■ 

We are interested now in the expression of the FT for 
the most general correlated gaussian 

h KLM {u,A;ir) = (n\ip KLM (u, A)) = 

(2n)- m / 2 [ e-*-"f KLU {u, A- x) dx = (66) 



(27T)" 



3N/2 / -xAx-iiv-xi- ■ \2K 



ux\ y LM (ux)dx. 



The principle for the proof is quite similar to what was 
proposed previously; the same changes of variables occur; 
the only difference concerns the integration over the Z\ 
variable that must be done now with help of the most 
general formula ([58]) . 

The FT for the general correlated gaussian takes the 
form 



h KLM (u,A;Tr) 



{-i) L K\ {uA^u) 



2£+3A72( det 4)3/2 



(67) 



xL 



K 



£+1/2 / (uA 1 TT) 2 



AuA- 



e 4 



n-A 1 n 



In this case, the FT of a general correlated gaus- 
sian is not proportional to a general correlated 
gaussian. 



However, using the series expansion for the Laguerre 
polynomial 



K{z) = £(-i)' 



T(n + a + I) 



(n — m)\T(m + a + 1) m! 



, (68) 



it follows that the FT of a general correlated gaus- 
sian is a linear combination of general correlated 
gaussians. Explicitly we have the property 



h KLM {u,A;-7r) 



( 2 3iv/2 det 4)3/2 K'\{K- K') 



(2K + 2L+1)!! (uA- 1 u) K - K ' , A' 1 

(2K' + 2L + 1)\\ 2^+K'+l 1k-lm{A u,—;tv). 

(69) 



VIII. RELATIVISTIC KINETIC ENERGY 

In many cases, the potential, or any operator, depends 
on the momenta 7r^ — —id/dxi. If this dependence 
occurs through an integer power of the momenta, the 
method explained above, for example in section HVi ap- 
plies. However this is not always the case. For such a 
delicate situation, there is only one remedy : to work in 
momentum representation. 

Let us suppose that the operator O depends on 
an argument of type wtt. Then the calculation of 
the matrix element (iPk'lm(u' , A')\0\i{}klm(u' , A')) = 
(hK>LM{u\A']Tv)\0(wTr)\h KLM (u,A;ir)) can be per- 
formed using the expansion (f6"9"|) . The dynamical element 
to be calculated has a form like 

-i A'^ 1 - A^ 1 
{fK'LKi{A' u, ;7r)|6(u)7r)|/K LM (^~ 1 w, —p ;tt)>. 

(70) 

This last element is then computed by application of Eq. 
(|47|) (or Eq. (f36|) ) with the trivial changes u — > A^u, 
A -> A-74, u' -> A'~V, A 1 -> A'~ X /A. 

As an example of application, we consider the relativis- 
tic kinetic energy operator 



Tn 



N+l 

E 



(71) 



The separation of the center of mass motion cannot be 
done properly in this case. The usual way to consider the 
intrinsic motion for the system is to work in the center 
of mass frame (which is well defined even in the spe- 
cial case where all particles have a vanishing mass); this 
means that we put P = in the subsequent formulae. 
In particular the momentum pi relative to the particle i 
writes (see relation ©) 



Pi 



N 

E^ 

j=i 



(72) 
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with 



W _ 



(i) 



= if j < i - 1, 
= -1 if j = * - 1, 
= li J > i- 



(73) 



In order to present simple formulae, let us focus on the 
special case K = K' = 0. Then 

(iI> lm{u',A')\T r \iIj 0L m{u,A)) = (74) 

N+l , 

E { h 0LM (u', A'; tt)| |l>«7r| 2 + mf \h 0LM (u, A; n) 



where K n {z) is the modified Bessel function. The integral 
for higher I can be obtained by the following recursion 
formula 



1 (21-1)0 + 1-2/3) 



All the above formulae allows to get a very efficient 
way to calculate the matrix element of the relativistic 
energy operator for the few-body problem expressed in 
terms of correlated gaussians. 

IX. CONCLUDING REMARKS 



The computation of the matrix element in the summation 
is performed first using Eq. (|6"5)) and then Eq. Ijl5]l with 
the appropriate parameters. 
The final result looks like 



&olm{u',A')\T r \iP olm (u,A)) 



\r N+l 

E z i 3/2 



-in 



E 



£j2»(2n + l)!(L-n)l 
x^(2n + 2, i L, mi ) 



(75) 



1 2Z i7i7 A 
P / 



L-r 



Again, the matrix element is proportional to the over- 
lap Ml- The p parameter is defined as in formula (fT8"|) 
while 



Z, = U®A!B- 1 AU®\ 



U^AB~ 1 u' 



U^A'B-h 



7i 



and the integral for the potential is given by 



F(k,A,m)= due- Au u k y / u 2 +' 



(76) 



(77) 



In fact, we need in Eq. (|75[) . the values of T(k,A,m) 
for even values of A; > 2. We calculated this function 
in the following way. Let us define for simplicity X; = 
jF(2Z + 2, A, m) and introduce the parameter (5 = Am 2 /2. 
The first values can be found analytically 



To = 
Ti = 



1 

2^3 



(78a) 



pe f3 [PK Q {p) + (2-P)Kim (78b) 



In this paper, we propose a number of new formulae 
concerning correlated gaussians that cannot be found in 
SV. Our formulation deals with general correlated gaus- 
sians (arbitrary number of particles, arbitrary orbital an- 
gular momentum, arbitrary value of the quantum num- 
ber K), but is limited to the natural parity states and 
to operators that do not mix spin and space degrees of 
freedom. 

We present new formulae relative to central poten- 
tials. Instead of integrals containing Hermitc polyno- 
mials, we proposed new formulae with simpler integrals 
which, most of the time, can be evaluated analytically. 
Moreover, we derived another formulation for the ma- 
trix elements which is more efficient than that presented 
in SV. All these formulae can be simplified a lot in the 
special case K = K 1 = 0. 

But the most interesting point of this paper is the 
derivation of the Fourier transform of correlated gaus- 
sians. In the peculiar case K = 0, we showed that 
the Fourier transform is just proportional to a correlated 
gaussian, but with renormalized parameters. In the gen- 
eral case, we showed that the Fourier transform is a com- 
binaison of general correlated gaussians. 

These new formulae allow the evaluation of the matrix 
elements of a relativistic kinetic energy operator. Since 
this kind of operator is more and more introduced in real- 
istic situations, the results of this paper are of primordial 
importance for future few-body calculations. 
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